{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "<p align=\"center\">\n",
    "    <img src=\"https://github.com/GeostatsGuy/GeostatsPy/blob/master/TCG_color_logo.png?raw=true\" width=\"220\" height=\"240\" />\n",
    "\n",
    "</p>\n",
    "\n",
    "### Monte Carlo Simulation in Python \n",
    "\n",
    "#### Michael Pyrcz, Professor, The University of Texas at Austin \n",
    "\n",
    "##### [Twitter](https://twitter.com/geostatsguy) | [GitHub](https://github.com/GeostatsGuy) | [Website](http://michaelpyrcz.com) | [GoogleScholar](https://scholar.google.com/citations?user=QVZ20eQAAAAJ&hl=en&oi=ao) | [Book](https://www.amazon.com/Geostatistical-Reservoir-Modeling-Michael-Pyrcz/dp/0199731446) | [YouTube](https://www.youtube.com/channel/UCLqEr-xV-ceHdXXXrTId5ig)  | [LinkedIn](https://www.linkedin.com/in/michael-pyrcz-61a648a1)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Monte Carlo Simulation\n",
    "\n",
    "Here's two dashboards to demonstration of Monte Carlo simulation for subsurface uncertainty modeling workflows. This should help you get started with:\n",
    "\n",
    "* building data science workflows integrate uncertainty sources, i.e, math with distributions\n",
    "* using more advanced data science workflows like bootstrap and Markov chain Monte Carlo (McMC)\n",
    "\n",
    "I have recorded a walk-through of this interactive dashboard in my [Data Science Interactive Python Demonstrations](https://www.youtube.com/playlist?list=PLG19vXLQHvSDy26fM3hDLg3VCU7U5BGZl) series on my [YouTube](https://www.youtube.com/@GeostatsGuyLectures) channel.\n",
    "\n",
    "* Join me for walk-through of this dashboard [09 Data Science Interactive: Monte Carlo Simulation](TBD). I'm stoked to guide you and share observations and things to try out!   \n",
    "\n",
    "* I have a lecture on [Monte Carlo Simulation](https://www.youtube.com/watch?v=Qb8TsSINpnU&list=PLG19vXLQHvSB-D4XKYieEku9GQMQyAzjJ&index=13&t=3s) as part of my [Data Analytics and Geostatistics](https://www.youtube.com/playlist?list=PLG19vXLQHvSB-D4XKYieEku9GQMQyAzjJ) course. Note, for all my recorded lecture the interactive and well-documented workflow demononstrations are available on my GitHub repository [GeostatsGuy's Python Numerical Demos](https://github.com/GeostatsGuy/PythonNumericalDemos).\n",
    "\n",
    "* Also, I have a lecture on [Markov chain Monte Carlo](https://www.youtube.com/watch?v=7QX-yVboLhk&list=PLG19vXLQHvSC2ZKFIkgVpI9fCjkN38kwf&index=34) as part of Bayesian regression in my [Machine Learning](https://www.youtube.com/playlist?list=PLG19vXLQHvSC2ZKFIkgVpI9fCjkN38kwf) course.\n",
    "\n",
    "#### Monte Carlo Simulation\n",
    "\n",
    "Definition: random sampling from a distribution.\n",
    "\n",
    "Procedure: \n",
    "\n",
    "1. Model the representative distribution (CDF), $F_X(x)$\n",
    "2. Draw a random value from a uniform [0,1] distribution (p-value), $p^\\ell$\n",
    "3. Apply the inverse of the CDF to calculate the associated realization, $F^{-1}_X(p^\\ell)$\n",
    "\n",
    "In practice, Monte Carlo Simulation (MCS) refers to the workflow with multiple realizations drawn to buld an uncertainty model. \n",
    "\n",
    "\\begin{equation}\n",
    "x^\\ell = F^{-1}_X(p^\\ell),  \\, \\forall \\, \\ell = 1,\\ldots, L\n",
    "\\end{equation}\n",
    "\n",
    "where $x^\\ell$ is the realization of the variable $X$ drawn from its CDF, $F_X(x)$, with cumulative probability, p-value, $p^\\ell$.  \n",
    "\n",
    "#### Monte Carlo Simulation for Uncertainty Modeling\n",
    "\n",
    "It would be trivial to apply MCS to a single random variable, after many realizations one would get back the original distribution. But this general method allows us to do math with distributions, random variables, to calculate uncertainty models! \n",
    "\n",
    "Definitions:\n",
    "\n",
    "* **random variable** - a feature's value at a time or location is unknown and take on a range of possible values, $X_1$, usually defined by a PDF, $f_{X_1}(x_1)$, or CDF, $F_{X_1}(x_1)$ \n",
    "\n",
    "* **data value** - a sample observation at a specific location or time for a feature, $x_1\\left(\\bf{u}_{\\alpha} \\right)$ \n",
    "\n",
    "* **realization** - a random sample drawn from a random variable, $x_{\\alpha}^{\\ell}$ \n",
    "\n",
    "With MCS we can solve problems like this:\n",
    "\n",
    "$Y = X_1 + X_2 + X_3$ given $X_1$, $X_2$, and $X_3$ are random variables with any distribution, parametric or empirical (a set of values known as a reference distribution).\n",
    "\n",
    "#### Monte Carlo Simulation for Uncertainty Modeling General Workflow\n",
    "\n",
    "The general approach is to:\n",
    "\n",
    "1. Model all distributions for the input, variables of interest $F_{X_1},\\ldots,F_{X_m}$.\n",
    "2. For each realization draw $p^\\ell_{1},\\ldots,p^\\ell_{m}$, p-values\n",
    "3. Apply the inverse of each distribution to calculate a realization of each variable, $x^\\ell_j = F_{X^\\ell_j}^{-1}(p^\\ell_j),  \\, \\forall \\, j = 1,\\ldots$, $m$ variables.\n",
    "4. Apply each set of variables for a $\\ell$ realization to the transfer function to calculate the ouptput realization, $y^\\ell = F(x_1^\\ell,\\ldots,x_m^\\ell)$.\n",
    "\n",
    "Monte Carlo Simulation (MCS) is extremely powerful\n",
    "\n",
    "* Possible to easily simulate uncertainty models for complicated systems \n",
    "* Simulations are conducted by drawing values at random from specified uncertainty distributions for each variable\n",
    "* A single realization of each variable, $X_1^\\ell, X_2^\\ell,\\ldots,X_m^\\ell$ is applied to the transfer function to calculate the realization of the variable of interest (output, decision criteria):\n",
    "\n",
    "\\begin{equation}\n",
    "Y^\\ell = f(X_1^\\ell,\\ldots,X_m^\\ell), \\, \\forall \\, \\ell = 1,\\ldots, L\n",
    "\\end{equation}\n",
    "\n",
    "* The MCS method builds empirical uncertainty models by random sampling\n",
    "\n",
    "#### How many realizations, $L$?\n",
    "\n",
    "The answer is enough realizations! \n",
    "\n",
    "* If the MCS computational cost is low then **many** is the right answer. \n",
    "\n",
    "* If too few realizations are calculated then the summary statistics and the entire CDF of the output, decision criteria may be incorrect. This is caused by fluctuations due to not enough samples (see the 'Law of Small Numbers'). \n",
    "\n",
    "#### The Simulation Perspective\n",
    "\n",
    "There are many uncertainty operations that we can work out analytically, for which we don't need MCS:\n",
    "\n",
    "* $X_1 + c$\n",
    "* $c \\cdot X_1$\n",
    "* $X_1 + X_2$, given $X_1$ and $X_2$ are Gaussian distributed\n",
    "\n",
    "But consider the following:\n",
    "\n",
    "* $X_1 \\times X_2$, given $X_1$ and $X_2$ could have any distribution\n",
    "* what is the probability of winning the game of solitaire?\n",
    "\n",
    "With MCS we simulate the system, sample $L$ outcomes and then summarize!\n",
    "\n",
    "#### Limitations\n",
    "\n",
    "The MCS method above assumes:\n",
    "1. **representativity** - the distributions are representative\n",
    "2. **independence** - the random variables are independent of eachother\n",
    "3. **stationarity** - all realizations for each variable are from the same distribution, e.g., no change over time nor space \n",
    "  \n",
    "#### Getting Started\n",
    "\n",
    "Here's the steps to get setup in Python with the GeostatsPy package:\n",
    "\n",
    "1. Install Anaconda 3 on your machine (https://www.anaconda.com/download/). \n",
    "\n",
    "That's all!\n",
    "\n",
    "#### Load the Required Libraries\n",
    "\n",
    "We will also need some standard Python packages. These should have been installed with Anaconda 3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "supress_warnings = True                                   # supress warnings?\n",
    "import numpy as np                                        # ndarrys for gridded data\n",
    "import matplotlib.pyplot as plt                           # for plotting\n",
    "from matplotlib.ticker import (MultipleLocator, AutoMinorLocator) # control of axes ticks\n",
    "from scipy import stats                                   # summary statistics\n",
    "from ipywidgets import interactive                        # widgets and interactivity\n",
    "from ipywidgets import widgets                            \n",
    "from ipywidgets import Layout\n",
    "from ipywidgets import Label\n",
    "from ipywidgets import VBox, HBox\n",
    "cmap = plt.cm.inferno                                     # default color bar, no bias and friendly for color vision defeciency\n",
    "plt.rc('axes', axisbelow=True)                            # grid behind plotting elements\n",
    "if supress_warnings == True:\n",
    "    import warnings                                       # supress any warnings for this demonstration\n",
    "    warnings.filterwarnings('ignore')                  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you get a package import error, you may have to first install some of these packages. This can usually be accomplished by opening up a command window on Windows and then typing 'python -m pip install [package-name]'. More assistance is available with the respective package docs.  \n",
    "\n",
    "#### Declare Functions\n",
    "\n",
    "I just added a convenience function for adding major and minor gridlines."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def add_grid():\n",
    "    plt.gca().grid(True, which='major',linewidth = 1.0); plt.gca().grid(True, which='minor',linewidth = 0.2) # add y grids\n",
    "    plt.gca().tick_params(which='major',length=7); plt.gca().tick_params(which='minor', length=4)\n",
    "    plt.gca().xaxis.set_minor_locator(AutoMinorLocator()); plt.gca().yaxis.set_minor_locator(AutoMinorLocator()) # turn on minor ticks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Simple Monte Carlo Simulation Demonstration\n",
    "\n",
    "Now let's set up our first dashboard."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "l = widgets.Text(value='                                          Monte Carlo Simulation Demonstration, Michael Pyrcz, Professor, The University of Texas at Austin',layout=Layout(width='950px', height='30px'))\n",
    "\n",
    "L = widgets.IntSlider(min=1, max = 40, value = 1, description = '$L$:',orientation='horizontal',layout=Layout(width='230px', height='50px'),continuous_update=False)\n",
    "L.style.handle_color = 'gray'\n",
    "\n",
    "dist1 = widgets.Dropdown(\n",
    "    options=['Uniform','Triangular','Gaussian'],\n",
    "    value='Gaussian',\n",
    "    description='$X_1$:',\n",
    "    disabled=False,\n",
    "    layout=Layout(width='200px', height='30px')\n",
    ")\n",
    "min1 = widgets.FloatSlider(min=0.0, max = 100.0, value = 10.0, description = 'Min',orientation='horizontal',layout=Layout(width='230px', height='50px'),continuous_update=False)\n",
    "min1.style.handle_color = 'gray'\n",
    "max1 = widgets.FloatSlider(min=0.0, max = 100.0, value = 90.0, description = 'Max',orientation='horizontal',layout=Layout(width='230px', height='50px'),continuous_update=False)\n",
    "max1.style.handle_color = 'gray'\n",
    "\n",
    "ui = widgets.HBox([L,dist1,min1,max1],kwargs = {'justify_content':'center'})\n",
    "ui1 = widgets.VBox([l,ui],)\n",
    "\n",
    "def make_sample_cdist(dist,zmin,zmax,L):\n",
    "    cdf = np.linspace(0.001,0.999,1000)\n",
    "    pvals = np.random.random(size = L)\n",
    "    if dist == 'Triangular':  \n",
    "        distribution = stats.triang(loc = zmin,c=0.5,scale = zmax-zmin)\n",
    "    if dist == 'Uniform':        \n",
    "        distribution = stats.uniform(loc = zmin,scale = zmax-zmin)\n",
    "    if dist == 'Gaussian':\n",
    "        mean = (zmax + zmin)*0.5; stdev = (zmax - zmin)/6.0\n",
    "        distribution = stats.norm(loc = mean,scale = stdev)\n",
    "    \n",
    "    cdfx = distribution.ppf(cdf)\n",
    "    sample = distribution.ppf(pvals)\n",
    "    \n",
    "    return sample, pvals, cdfx, cdf \n",
    "        \n",
    "def f_make1(L,dist1,min1,max1): \n",
    "    np.random.seed(seed = 73073)\n",
    "    sample, pvals, cdfx, cdf = make_sample_cdist(dist1,min1,max1,L)\n",
    "    \n",
    "    plt.subplot(121)\n",
    "    plt.plot(cdfx,cdf,'--',color='black',linewidth = 3,zorder=10); add_grid()\n",
    "    plt.xlim(0,100); plt.ylim(0,1.0); plt.xlabel(\"$X_1$\"); plt.title(\"Cumulative Distribution Function, $X_1$\"); plt.ylabel('Cumulative Probability')\n",
    "     \n",
    "    for l in range(0,L):\n",
    "        alpha = max(0.02,(6 - (L-l) )/5)\n",
    "        dhead = 1;lhead = 0.02        \n",
    "        if l == L-1:\n",
    "            plt.plot([0.0,sample[l]],[pvals[l],pvals[l]],color='red',alpha = alpha,lw=1,zorder=1)\n",
    "            plt.plot([sample[l],sample[l]],[pvals[l],lhead+0.02],color='red',alpha = alpha,lw=1,zorder=1)\n",
    "            head = plt.Polygon([[sample[l],0.02],[sample[l]-dhead,lhead+0.02],[sample[l]+dhead,lhead+0.02],[sample[l],0.02]], color='red',alpha = alpha,zorder=1)\n",
    "            plt.annotate(int(sample[l]),[sample[l]+1,0.07],rotation=90.0,color='red',alpha = alpha)\n",
    "            plt.scatter(sample[l],0.01,s=30,color='red',edgecolor='black',alpha=1.0)\n",
    "            plt.gca().add_patch(head)\n",
    "        else:\n",
    "            plt.plot([0.0,sample[l]],[pvals[l],pvals[l]],color='darkorange',alpha = alpha,lw=1,zorder=1)\n",
    "            plt.plot([sample[l],sample[l]],[pvals[l],lhead+0.02],color='darkorange',alpha = alpha,lw=1,zorder=1)\n",
    "            head = plt.Polygon([[sample[l],0.02],[sample[l]-dhead,lhead+0.02],[sample[l]+dhead,lhead+0.02],[sample[l],0.02]], color='darkorange',alpha = alpha,zorder=1)\n",
    "            plt.annotate(int(sample[l]),[sample[l]+1,0.07],rotation=90.0,color='darkorange',alpha = alpha)\n",
    "            plt.scatter(sample[l],0.01,s=10,color='yellow',edgecolor='black',alpha=1.0)\n",
    "            plt.gca().add_patch(head)            \n",
    "\n",
    "    plt.subplot(122)\n",
    "    plt.hist(sample,density= False,bins=np.linspace(0,100,20),weights=None,color='red',alpha=1.0,edgecolor='black',zorder=1)\n",
    "    if l > 0:\n",
    "        plt.hist(sample[:-1],density= False,bins=np.linspace(0,100,20),weights=None,color='darkorange',alpha=0.8,edgecolor='black',zorder=5)\n",
    "    \n",
    "    for l in range(0,L):\n",
    "        if l == L-1:\n",
    "            plt.scatter(sample[l],0.1,s=30,color='red',edgecolor='black',alpha=1.0,zorder=l+10)\n",
    "        else:\n",
    "            plt.scatter(sample[l],0.1,s=10,color='yellow',edgecolor='black',alpha=1.0,zorder=l+10)\n",
    "    \n",
    "    plt.xlabel(r\"$x_{\\alpha}, \\alpha = 1,\\ldots,L$\"); plt.title(r\"Monte Carlo Simulation Realizations, $x_{\\alpha}, \\alpha = 1,\\ldots,L$\"); plt.ylabel('Frequency')\n",
    "    plt.xlim([0,100]); plt.ylim([0,10]); add_grid()\n",
    "    \n",
    "    plt.subplots_adjust(left=0.0, bottom=0.0, right=2.0, top=1.2, wspace=0.3, hspace=0.2)\n",
    "    plt.show()    \n",
    "\n",
    "interactive_plot1 = widgets.interactive_output(f_make1, {'L':L,'dist1':dist1,'min1':min1,'max1':max1})\n",
    "interactive_plot1.clear_output(wait = True)                # reduce flickering by delaying plot updating  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "495be71ffccd45098a8f56d8d90d5288",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(Text(value='                                          Monte Carlo Simulation Demonstration, Mic…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "bd669e6ce2154b04a1e90dc90c1c365b",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display(ui1, interactive_plot1)                            # display the interactive plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Monte Carlo Simulation Demonstration\n",
    "\n",
    "Now let's set up our second dashboard for adding or multiplying distributions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# interactive calculation of the random sample set (control of source parametric distribution and number of samples)\n",
    "l = widgets.Text(value='                                          Monte Carlo Simulation Demonstration, Michael Pyrcz, Professor, The University of Texas at Austin',layout=Layout(width='950px', height='30px'))\n",
    "\n",
    "operator = widgets.RadioButtons(options=['Add','Mult'],description='Operator:',disabled=False,layout=Layout(width='230px', height='100px'))\n",
    "\n",
    "L = widgets.IntSlider(min=1, max = 10000, value = 1, description = '$L$:',orientation='horizontal',layout=Layout(width='230px', height='50px'),continuous_update=False)\n",
    "L.style.handle_color = 'gray'\n",
    "\n",
    "uiL = widgets.VBox([L,operator])\n",
    "\n",
    "dist1 = widgets.Dropdown(\n",
    "    options=['Uniform','Triangular','Gaussian'],\n",
    "    value='Gaussian',\n",
    "    description='$X_1$:',\n",
    "    disabled=False,\n",
    "    layout=Layout(width='200px', height='30px')\n",
    ")\n",
    "min1 = widgets.FloatSlider(min=0.0, max = 100.0, value = 10.0, description = 'Min',orientation='horizontal',layout=Layout(width='230px', height='50px'),continuous_update=False)\n",
    "min1.style.handle_color = 'blue'\n",
    "max1 = widgets.FloatSlider(min=0.0, max = 100.0, value = 30.0, description = 'Max',orientation='horizontal',layout=Layout(width='230px', height='50px'),continuous_update=False)\n",
    "max1.style.handle_color = 'blue'\n",
    "ui1 = widgets.VBox([dist1,min1,max1],kwargs = {'justify_content':'center'}) \n",
    "\n",
    "dist2 = widgets.Dropdown(\n",
    "    options=['Triangular', 'Uniform', 'Gaussian'],\n",
    "    value='Gaussian',\n",
    "    description='$X_2$:',\n",
    "    disabled=False,\n",
    "    layout=Layout(width='200px', height='30px')\n",
    ")\n",
    "min2 = widgets.FloatSlider(min=0.0, max = 100.0, value = 10.0, description = 'Min',orientation='horizontal',layout=Layout(width='230px', height='50px'),continuous_update=False)\n",
    "min2.style.handle_color = 'red'\n",
    "max2 = widgets.FloatSlider(min=0.0, max = 100.0, value = 30.0, description = 'Max',orientation='horizontal',layout=Layout(width='230px', height='50px'),continuous_update=False)\n",
    "max2.style.handle_color = 'red'\n",
    "ui2 = widgets.VBox([dist2,min2,max2],kwargs = {'justify_content':'center'})\n",
    "\n",
    "dist3 = widgets.Dropdown(\n",
    "    options=['Triangular', 'Uniform', 'Gaussian'],\n",
    "    value='Gaussian',\n",
    "    description='$X_3$:',\n",
    "    disabled=False,\n",
    "    layout=Layout(width='200px', height='30px')\n",
    ")\n",
    "min3 = widgets.FloatSlider(min=0.0, max = 100.0, value = 10.0, description = 'Min',orientation='horizontal',layout=Layout(width='230px', height='50px'),continuous_update=False)\n",
    "min3.style.handle_color = 'yellow'\n",
    "max3 = widgets.FloatSlider(min=0.0, max = 100.0, value = 30.0, description = 'Max',orientation='horizontal',layout=Layout(width='230px', height='50px'),continuous_update=False)\n",
    "max3.style.handle_color = 'yellow'\n",
    "ui3 = widgets.VBox([dist3,min3,max3],kwargs = {'justify_content':'center'})\n",
    "\n",
    "ui = widgets.HBox([uiL,ui1,ui2,ui3])\n",
    "ui2 = widgets.VBox([l,ui],)\n",
    "\n",
    "def make_dist(dist,zmin,zmax,L):\n",
    "    \n",
    "    if dist == 'Triangular':\n",
    "        z = np.random.triangular(left=zmin, mode=(zmax+zmin)*0.5, right=zmax, size=L)\n",
    "        pdf = stats.triang.pdf(np.linspace(0.0,100.0,1000), loc = zmin, c = 0.5, scale = zmax-zmin)* 2 * L \n",
    "    if dist == 'Uniform':\n",
    "        z = np.random.uniform(low=zmin, high=zmax, size=L)\n",
    "        pdf = stats.uniform.pdf(np.linspace(0.0,100.0,1000), loc = zmin, scale = zmax-zmin) * 2 * L\n",
    "    if dist == 'Gaussian':\n",
    "        mean = (zmax + zmin)*0.5; sd = (zmax - zmin)/6.0\n",
    "        z = np.random.normal(loc = mean, scale = sd, size=L)\n",
    "        pdf = stats.norm.pdf(np.linspace(0.0,100.0,1000), loc = mean, scale = sd) * 2 * L\n",
    "    return z, pdf\n",
    "        \n",
    "def f_make(L,operator,dist1,min1,max1,dist2,min2,max2,dist3,min3,max3): \n",
    "    np.random.seed(seed = 73073)\n",
    "    x1, pdf1 = make_dist(dist1,min1,max1,L)\n",
    "    x2, pdf2 = make_dist(dist2,min2,max2,L)\n",
    "    x3, pdf3 = make_dist(dist3,min3,max3,L)\n",
    "\n",
    "    xvals = np.linspace(0.0,100.0,1000)\n",
    "    plt.subplot(241)\n",
    "    plt.hist(x1,density = False,bins=np.linspace(0,100,50),weights=None,color='blue',alpha=0.7,edgecolor='grey')\n",
    "    plt.plot(xvals,pdf1,'--',color='black',linewidth = 3)\n",
    "    plt.xlim(0,100); plt.xlabel(\"$X_1$\"); plt.title(\"First Predictor Feature, $X_1$\"); plt.ylabel('Frequency')\n",
    " \n",
    "    plt.subplot(242)\n",
    "    plt.hist(x2,density = False,bins=np.linspace(0,100,50),weights=None,color='red',alpha=0.7,edgecolor='grey')\n",
    "    plt.plot(xvals,pdf2,'--',color='black',linewidth = 3)\n",
    "    plt.xlim(0,100); plt.xlabel(\"$X_1$\"); plt.title(\"Second Predictor Feature, $X_2$\"); plt.ylabel('Frequency')\n",
    " \n",
    "    plt.subplot(243)\n",
    "    plt.hist(x3,density = False,bins=np.linspace(0,100,50),weights=None,color='yellow',alpha=0.7,edgecolor='grey')\n",
    "    plt.plot(xvals,pdf3,'--',color='black',linewidth = 3)\n",
    "    plt.xlim(0,100); plt.xlabel(\"$X_1$\"); plt.title(\"Third Predictor Feature, $X_3$\"); plt.ylabel('Frequency')\n",
    " \n",
    "    y = np.zeros(L)\n",
    "    ymin = 0.0\n",
    "    if operator == \"Add\":\n",
    "        y = x1 + x2 + x3\n",
    "    elif operator == \"Mult\":\n",
    "        y = x1 * x2 * x3\n",
    "        \n",
    "    ymax = max(round((np.max(y)+50)/100)*100,100) # round up to nearest hundreds to avoid the chart jumping around\n",
    "    \n",
    "    plt.subplot(244)\n",
    "    plt.hist(y,density = False,bins=np.linspace(ymin,ymax,50),weights=None,color='black',alpha=0.5,edgecolor='black')\n",
    "    plt.xlabel(\"$Y$\"); plt.ylabel('Frequency'); plt.xlim(ymin,ymax)\n",
    "    if operator == \"Add\":\n",
    "        plt.title(\"Response Feature, $y = X_1 + X_2 + X_3$\")\n",
    "    else:\n",
    "        plt.title(r\"Response Feature, $y = X_1 \\times X_2 \\times X_3$\")\n",
    "    \n",
    "    plt.subplots_adjust(left=0.0, bottom=0.0, right=2.0, top=1.2, wspace=0.3, hspace=0.2)\n",
    "    plt.show()    \n",
    "\n",
    "interactive_plot = widgets.interactive_output(f_make, {'L':L,'operator':operator,'dist1':dist1,'min1':min1,'max1':max1,'dist2':dist2,'min2':min2,'max2':max2,'dist3':dist3,'min3':min3,'max3':max3})\n",
    "interactive_plot.clear_output(wait = True)                # reduce flickering by delaying plot updating    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Monte Carlo Simulation Demonstration\n",
    "\n",
    "* specify the distributions for 3 Random Variables, $X_1$, $X_2$, and $X_3$ and select the operator $Y = f(X_1,X_2,X_3)$\n",
    "\n",
    "* observe the distribution of the resulting Monte Carlos Simulation realization histograms of $x_1^{\\ell}$, $x_2^{\\ell}$, $x_3^{\\ell}$, and $y^{\\ell}$\n",
    "\n",
    "#### Michael Pyrcz, Professor, The University of Texas at Austin \n",
    "\n",
    "##### [Twitter](https://twitter.com/geostatsguy) | [GitHub](https://github.com/GeostatsGuy) | [Website](http://michaelpyrcz.com) | [GoogleScholar](https://scholar.google.com/citations?user=QVZ20eQAAAAJ&hl=en&oi=ao) | [Book](https://www.amazon.com/Geostatistical-Reservoir-Modeling-Michael-Pyrcz/dp/0199731446) | [YouTube](https://www.youtube.com/channel/UCLqEr-xV-ceHdXXXrTId5ig)  | [LinkedIn](https://www.linkedin.com/in/michael-pyrcz-61a648a1) | [GeostatsPy](https://github.com/GeostatsGuy/GeostatsPy)\n",
    "\n",
    "### The Inputs\n",
    "\n",
    "* **$L$**: number of realizations, **Operator**: addition for $Y = X_1 + X_2 + X_3$, multiplication for $Y = X_1 \\times X_2 \\times X_3$\n",
    "\n",
    "* **$X_1$, $X_2$, and $X_3$**: distribution type, min and max. Assume mode or mean is centered and 3 st.dev. for Gaussian"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "a1940971104b4fc29ae97f9f0d53c50a",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(Text(value='                                          Monte Carlo Simulation Demonstration, Mic…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "c9901773756f4d0abc847746f480a73e",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output(outputs=({'output_type': 'display_data', 'data': {'text/plain': '<Figure size 432x288 with 4 Axes>', 'i…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display(ui2, interactive_plot)                            # display the interactive plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Comments\n",
    "\n",
    "This was a basic demonstration of Monte Carlo simulation for uncertainty analysis. A lot more could be done, for example, more complicated transfer functions and a combination of non-parametric and parametric distributions. Also, one could integrate relationships between the variables (we assumed independent here).\n",
    "\n",
    "Note, I have other demonstrations with McMC methods, also!\n",
    "  \n",
    "I hope this was helpful,\n",
    "\n",
    "*Michael*\n",
    "\n",
    "#### The Author:\n",
    "\n",
    "### Michael Pyrcz, Professor, The University of Texas at Austin \n",
    "*Novel Data Analytics, Geostatistics and Machine Learning Subsurface Solutions*\n",
    "\n",
    "With over 17 years of experience in subsurface consulting, research and development, Michael has returned to academia driven by his passion for teaching and enthusiasm for enhancing engineers' and geoscientists' impact in subsurface resource development. \n",
    "\n",
    "For more about Michael check out these links:\n",
    "\n",
    "#### [Twitter](https://twitter.com/geostatsguy) | [GitHub](https://github.com/GeostatsGuy) | [Website](http://michaelpyrcz.com) | [GoogleScholar](https://scholar.google.com/citations?user=QVZ20eQAAAAJ&hl=en&oi=ao) | [Book](https://www.amazon.com/Geostatistical-Reservoir-Modeling-Michael-Pyrcz/dp/0199731446) | [YouTube](https://www.youtube.com/channel/UCLqEr-xV-ceHdXXXrTId5ig)  | [LinkedIn](https://www.linkedin.com/in/michael-pyrcz-61a648a1)\n",
    "\n",
    "#### Want to Work Together?\n",
    "\n",
    "I hope this content is helpful to those that want to learn more about subsurface modeling, data analytics and machine learning. Students and working professionals are welcome to participate.\n",
    "\n",
    "* Want to invite me to visit your company for training, mentoring, project review, workflow design and / or consulting? I'd be happy to drop by and work with you! \n",
    "\n",
    "* Interested in partnering, supporting my graduate student research or my Subsurface Data Analytics and Machine Learning consortium (co-PIs including Profs. Foster, Torres-Verdin and van Oort)? My research combines data analytics, stochastic modeling and machine learning theory with practice to develop novel methods and workflows to add value. We are solving challenging subsurface problems!\n",
    "\n",
    "* I can be reached at mpyrcz@austin.utexas.edu.\n",
    "\n",
    "I'm always happy to discuss,\n",
    "\n",
    "*Michael*\n",
    "\n",
    "Michael Pyrcz, Ph.D., P.Eng. Professor, Cockrell School of Engineering and The Jackson School of Geosciences, The University of Texas at Austin\n",
    "\n",
    "#### More Resources Available at: [Twitter](https://twitter.com/geostatsguy) | [GitHub](https://github.com/GeostatsGuy) | [Website](http://michaelpyrcz.com) | [GoogleScholar](https://scholar.google.com/citations?user=QVZ20eQAAAAJ&hl=en&oi=ao) | [Book](https://www.amazon.com/Geostatistical-Reservoir-Modeling-Michael-Pyrcz/dp/0199731446) | [YouTube](https://www.youtube.com/channel/UCLqEr-xV-ceHdXXXrTId5ig)  | [LinkedIn](https://www.linkedin.com/in/michael-pyrcz-61a648a1)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
